home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / nrpas13.zip / LFIT.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  81 lines
  1. PROCEDURE lfit(x,y,sig: glndata; ndata: integer; VAR a: glmma; mma: integer;
  2.       lista: gllista; mfit: integer; VAR covar: glcovar;
  3.       ncvm: integer; VAR chisq: real);
  4. (* Programs using routine LFIT must define the types
  5. TYPE
  6.    glndata = ARRAY [1..ndata] OF real;
  7.    glmma = ARRAY [1..mma] OF real;
  8.    gllista = ARRAY [1..mma] OF integer;
  9.    glcovar = ARRAY [1..ncvm,1..ncvm] OF real;
  10.    glnpbymp = ARRAY [1..ncvm,1..1] OF real;
  11. in the main routine. *)
  12. VAR
  13.    k,kk,j,ihit,i: integer;
  14.    ym,wt,sum,sig2i: real;
  15.    beta: glnpbymp;
  16.    afunc: glmma;
  17. BEGIN
  18.    kk := mfit+1;
  19.    FOR j := 1 TO mma DO BEGIN
  20.       ihit := 0;
  21.       FOR k := 1 TO mfit DO BEGIN
  22.          IF (lista[k] = j) THEN ihit := ihit+1
  23.       END;
  24.       IF (ihit = 0) THEN BEGIN
  25.          lista[kk] := j;
  26.          kk := kk+1
  27.       END ELSE IF (ihit > 1) THEN BEGIN
  28.          writeln('pause in routine LFIT');
  29.          writeln('improper permutation in LISTA'); readln
  30.       END
  31.    END;
  32.    IF (kk <> (mma+1)) THEN BEGIN
  33.       writeln('pause in routine LFIT');
  34.       writeln('improper permutation in LISTA'); readln
  35.    END;
  36.    FOR j := 1 TO mfit DO BEGIN
  37.       FOR k := 1 TO mfit DO BEGIN
  38.          covar[j,k] := 0.0
  39.       END;
  40.       beta[j,1] := 0.0
  41.    END;
  42.    FOR i := 1 TO ndata DO BEGIN
  43.       funcs(x[i],afunc,mma);
  44.       ym := y[i];
  45.       IF (mfit < mma) THEN BEGIN
  46.          FOR j := (mfit+1) TO mma DO BEGIN
  47.             ym := ym-a[lista[j]]*afunc[lista[j]]
  48.          END
  49.       END;
  50.       sig2i := 1.0/sqr(sig[i]);
  51.       FOR j := 1 TO mfit DO BEGIN
  52.          wt := afunc[lista[j]]*sig2i;
  53.          FOR k := 1 TO j DO BEGIN
  54.             covar[j,k] := covar[j,k]+wt*afunc[lista[k]]
  55.          END;
  56.          beta[j,1] := beta[j,1]+ym*wt
  57.       END
  58.    END;
  59.    IF (mfit > 1) THEN BEGIN
  60.       FOR j := 2 TO mfit DO BEGIN
  61.          FOR k := 1 TO j-1 DO BEGIN
  62.             covar[k,j] := covar[j,k]
  63.          END
  64.       END
  65.    END;
  66.    gaussj(covar,mfit,ncvm,beta,1,1);
  67.    FOR j := 1 TO mfit DO BEGIN
  68.       a[lista[j]] := beta[j,1]
  69.    END;
  70.    chisq := 0.0;
  71.    FOR i := 1 TO ndata DO BEGIN
  72.       funcs(x[i],afunc,mma);
  73.       sum := 0.0;
  74.       FOR j := 1 TO mma DO BEGIN
  75.          sum := sum+a[j]*afunc[j]
  76.       END;
  77.       chisq := chisq+sqr((y[i]-sum)/sig[i])
  78.    END;
  79.    covsrt(covar,ncvm,mma,lista,mfit)
  80. END;
  81.